﻿#pragma once

#include <algorithm>
#include <array>
#include <cmath>
#include <limits>
#include <numeric>
#include <tuple>
#include <type_traits>

namespace xy
{

namespace math
{

/**
 * @brief 数值常量类
 *
 */
struct numeric
{
    static constexpr double epsilon = 1e-12;
    static constexpr double distance = 1e-8;
    static constexpr double half = 1e-6;
    static constexpr std::size_t iter = 100;
    static constexpr double pi = 3.14159265358979323846;
    static constexpr double e = 2.71828182845904523536;
};

/**
 * @brief Hash 类模板，用于计算 std::tuple<T, T> 的哈希值
 *
 * @tparam T
 */
template <typename T>
    requires std::is_integral_v<T>
struct hash_tuple
{
    std::size_t operator()(const std::tuple<T, T> &key) const
    {
        auto &[x, y] = key;
        return (x * 2718281828) ^ (y * 3141592653);
    }
};

/**
 * @brief 标准比较器
 *
 */
struct compare
{
    using has_abs = decltype([](auto &&x) -> decltype(std::abs(x)) {});
    using has_distance = decltype([](auto &&x, auto &&y) -> decltype(x.distance(y)) {});

    template <typename T, typename U>
    constexpr auto less(const U &x, const U &y) const -> decltype(std::declval<has_abs>()(std::declval<U>()), bool{})
    {
        return std::abs(x - y) < std::numeric_limits<T>::epsilon();
    }

    template <typename T, typename U>
    constexpr auto less(const U &x, const U &y) const
        -> decltype(std::declval<has_distance>()(std::declval<U>(), std::declval<U>()), bool{})
    {
        return x.distance(y) < std::numeric_limits<T>::epsilon();
    }
};

/**
 * @brief 符号函数
 *
 * @tparam T
 * @param[in] t
 * @return int
 */
template <class T> inline int sign(T t)
{
    return t > 0 ? 1 : -1;
}

/**
 * @brief 阶乘函数（这里有一个巨大的坑，-1 * size_t 会得到负的 size_t，导致计算溢出，因此返回 int 类型而不是 size_t
 * 类型）
 *
 * @param[in] n
 * @return
 */
inline int factorial(std::size_t n)
{
    std::size_t i = 1;
    for (std::size_t k = 0; k < n; k++)
        i *= k + 1;
    return i;
}

/**
 * @brief 幂函数
 *
 * @tparam T
 * @param[in] base
 * @param[in] exp
 * @return requires
 */
template <typename T>
    requires std::is_integral_v<T>
inline T power(T base, T exp)
{
    T result = 1;
    for (std::size_t i = 0; i < exp; i++)
        result *= base;
    return result;
}

/**
 * @brief 排列数（这里有一个巨大的坑，-1 * size_t 会得到负的 size_t，导致计算溢出，因此返回 int 类型而不是 size_t 类型）
 *
 * @param[in] n
 * @param[in] i
 * @return
 */
inline int permutation(std::size_t n, std::size_t i)
{
    int j = 1;
    for (int k = n; k > n - i + 1; k--)
        j *= k;
    return j;
}

/**
 * @brief 组合数（这里有一个巨大的坑，-1 * size_t 会得到负的 size_t，导致计算溢出，因此返回 int 类型而不是 size_t 类型）
 * binomial(n, i) = binomial(n - 1, i) + binomial(n - 1, i - 1)
 * binomial(n, i) = binomial(n - 1, i - 1) * (n - i + 1) / i
 * 使用第二个公式防止整型溢出
 *
 * @param[in] n
 * @param[in] i
 * @return
 */
inline int binomial(std::size_t n, std::size_t i)
{
    if (i > n - i)
        i = n - i;

    int j = 1;
    for (std::size_t k = 1; k <= i; k++)
        j = j * (n - k + 1) / k;
    return j;
}

/**
 * @brief 裁剪
 *
 * @tparam T
 * @param t
 * @param min
 * @param max
 * @return T
 */
template <typename T> inline T clamp(T t, T min, T max)
{
    return std::min(std::max(t, min), max);
}

/**
 * @brief 角度转弧度
 *
 * @tparam T
 * @param[in] degrees
 * @return T
 */
template <typename T> inline T degrees_to_radians(T degrees)
{
    return degrees / 180.0 * numeric::pi;
}

/**
 * @brief 弧度转角度
 *
 * @tparam T
 * @param[in] radians
 * @return T
 */
template <typename T> inline T radians_to_degrees(T radians)
{
    return radians / numeric::pi * 180.0;
}

/**
 * @brief 基本数据内积
 *
 * @param[in] lhs
 * @param[in] rhs
 * @return double
 */
inline double dot(double lhs, double rhs)
{
    return lhs * rhs;
}

/**
 * @brief 求解 2 阶方程组
 *
 * @param[in] A
 * @param[in] b
 * @return true
 * @return false
 */
inline bool solve(double A[], double b[])
{
    double det = A[0] * A[3] - A[1] * A[2];
    if (std::abs(det) < std::numeric_limits<double>::epsilon())
        return false;

    auto b0 = (A[3] * b[0] - A[1] * b[1]) / det;
    auto b1 = (-A[2] * b[0] + A[0] * b[1]) / det;
    b[0] = b0;
    b[1] = b1;
    return true;
}

/**
 * @brief 计算 householder 变换
 *
 * @param[in] x
 * @param[in] n
 * @return double
 */
inline double householder(double x[], int n)
{
    // 使用无穷范数规范化防止溢出
    double eta = *std::max_element(x, x + n, [](double a, double b) { return std::abs(a) < std::abs(b); });
    std::transform(x, x + n, x, [eta](double a) { return a / eta; });
    double sigma = std::inner_product(x + 1, x + n, x + 1, 0.0);

    double beta = 0;
    if (std::abs(sigma) > std::numeric_limits<double>::epsilon())
    {
        double alpha = std::sqrt(x[0] * x[0] + sigma);
        if (x[0] <= 0)
            x[0] -= alpha;
        // 如果是正的分量，则用特殊算法减小舍入误差
        else
            x[0] = -sigma / (x[0] + alpha);

        // 对 beta 乘 v[0] * v[0] ，抵消规格化为 1 的影响
        beta = 2 * x[0] * x[0] / (sigma + x[0] * x[0]);
        std::transform(x, x + n, x, [gamma = x[0]](double a) { return a / gamma; });
    }
    return beta;
}

/**
 * @brief Givens 变换
 *
 * @param[in] v
 * @param[in] i
 * @param[in] j
 * @return std::tuple<double, double>
 */
inline std::tuple<double, double> givens(double v[], int i, int j)
{
    double sqrtX = std::sqrt(v[i] * v[i] + v[j] * v[j]);
    if (sqrtX == 0)
        return {1, 0};
    return {v[i] / sqrtX, v[j] / sqrtX};
}

} // namespace math

} // namespace xy
